function a = sa(N, NN, t)

  function res = Sa(X)
    % 2 x Sa(x)
    res = zeros(length(X), 1);
    for i=1:length(X)
      x = X(i);
      if (t == 2)
        if (x < 0.001)
          res(i) = (4/3).*x.^2+(-3/2).*x.^3+(5/3).*x.^5+(-12/7).*x.^6+(16/9).*x.^8+(-9/5).*x.^9+(11/6).*x.^11+(-24/13).*x.^12+(28/15).*x.^14+(-15/8).*x.^15+(17/9).*x.^17+(-36/19).*x.^18+(40/21).*x.^20;
        else
          res(i) = (1/3).*x.^(-1).*(3.^(1/2).*pi+(-6).*x.*(1+x).*(1+x+x.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*x.^(-1).*(2+x))+3.*log(1+x+x.^2));
        end
      elseif (t == 4)
        if (x < 0.001)
          res(i) = (-1/2)+(-1/3).*x.^2+(3/8).*x.^3+(-5/12).*x.^5+(3/7).*x.^6+(-4/9).*x.^8+(9/20).*x.^9+(-11/24).*x.^11+(6/13).*x.^12+(-7/15).*x.^14+(15/32).*x.^15+(-17/36).*x.^17+(9/19).*x.^18+(-10/21).*x.^20;
        else
          res(i) = (-1/12).*x.^(-1).*(3.^(1/2).*pi+6.*x.^3.*(1+x+x.^2).^(-1)+(-2).*3.^(1/2).*atan(3.^(-1/2).*x.^(-1).*(2+x))+3.*log(1+x+x.^2));
        end
      elseif (t == 5)
        if (x < 0.001)
          res(i) = (1/2)+x.^2+(-9/8).*x.^3+(5/4).*x.^5+(-9/7).*x.^6+(4/3).*x.^8+(-27/20).*x.^9+(11/8).*x.^11+(-18/13).*x.^12+(7/5).*x.^14+(-45/32).*x.^15+(17/12).*x.^17+(-27/19).*x.^18+(10/7).*x.^20;
        else
          res(i) = (1/2).*(((-2)+(-2).*x+x.^2).*(1+x+x.^2).^(-1)+(1/2).*x.^(-1).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*x.^(-1).*(2+x))+3.*log(1+x+x.^2)));
        end
      end
    end
  end
  z1 = chebpoints(NN, 0, 1)';
  D = chebop(NN, 0, 1);

  L = (diag(2+z1.^3) + diag(z1.*(z1.^3-1))*D)*D;

  b = 2.*z1.*Sa(z1);
  L(NN, :) = fliplr(eye(1,NN));
  b(NN) = 0;
  a1 = L\b;

  z1lin = ([0:NN-1]./(NN-1))';
  a1lin = chebinterp(a1, z1lin, 0, 1);
  a1lin = a1lin - a1lin(end);

% *****

  z2 = chebpoints(NN, 1, 2)';
  D = chebop(NN, 1, 2);

  L = (diag(2+z2.^3) + diag(z2.*(z2.^3-1))*D)*D;
  b = 2.*z2.*Sa(z2);
  L(1, :) = fliplr(eye(1,NN));
  b(1) = 0;
  a2 = L\b;

  z2lin = ([0:NN-1]./(NN-1)+1)';
  a2lin = chebinterp(a2, z2lin, 1, 2);
  a2lin = a2lin - a2lin(1);


  alin = [a1lin; a2lin(2:NN)];
  zlin = [z1lin; z2lin(2:NN)];
  alin = alin - alin(1);
  
  z = chebpoints(N, 0, 2);
  a = interp1(zlin, alin, z);

  plot(z, a)

  global res1;
  res1 = a;



end
